% This is the functions file of helicopter
% The full description of question3 is in question3.m
% structure --properties : parameters of an helicopter
%           --methods : functions for thrust,Cp,power,climb_rate,theta
% clear;clc
% Using object-oriented programming method to easily add the methods or
% properties and control the source code


classdef helicopter < handle
    properties(Constant)
        % environment
        g = 9.81; % gravity
        rho = 1.2256; % Air density at sea level
        
        % Rotor
        R = 12.8/2; % Main rotor radius 6.4m
        chrod =0.4; % Main rotor chord length(c)
        m_blade = 59.76; % Main rotor blade mass
        R_cq = 3.2; % Main rotor blade cg location
        R_t = 2.36/2; % Tail rotor diameter
        chrod_t = 0.18; % Tail rotor chord
        a = 5.5; % Main and tail rotors aerofoil lift curve slope
        theta1 = -9; % Main rotor blade washout
        omega = 360*pi/30; % Constant main rotor speed
        b = 4; % number of the blades        
        
        % Helicopter
        L = 13.33; % Fuselage length
        S_fp = 1.6; % Fuselage vertical flight equivalent flat plate drag area
        f = 1.8; % Fuselage forward flight equivalent flat plate drag area
        l_t = 7.66; % Distance between main rotor and tail rotor

        % Weights
        m_maxTakeoff = 5330; % Max take off
        m_maxLoad = 1360; % Max external load
        m_empty = 3291; % Empty weight
        m_maxFuel = 786; % Max fuel
        m_fuelLoad = 756; % individual fuel load
        m_payLoad = 910; % individual payload

        % Performance
        k1 = 1.15;
        k2 = 4.6;
        lamda_c=0; % initial value for climb inflow ratio

        % Engine Informance
        consu = 7.85413e-7 % Specific Fuel Consumption

    end
    methods(Static=true)
        function y = thrust(mu)
            % This function is aimed to calculate the thrust of a main
            % rotor
            g = helicopter.g;
            R = helicopter.R;
            S_fp = helicopter.S_fp;
            m_empty = helicopter.m_empty; 
            m_fuelLoad = helicopter.m_fuelLoad;
            m_payLoad = helicopter.m_payLoad;        

            m = m_empty + m_fuelLoad + m_payLoad;
            A = pi * R^2;
            W = m * g;
            T = W/(1-S_fp/A); % rotor thrust which accounts for the vertical drag
            y = T;
        end
    end
    methods (Static = true)
        % This function is aimed to calculate total helicopter Power coefficient
        function y = Cp(mu)
            g = helicopter.g;
            rho = helicopter.rho; 
            R = helicopter.R;
            chrod =helicopter.chrod; 
            m_blade = helicopter.m_blade; 
            R_cq = helicopter.R_cq; 
            R_t = helicopter.R_t; 
            chrod_t = helicopter.chrod_t; 
            a = helicopter.a; 
            theta1 = helicopter.theta1;
            omega = helicopter.omega;
            b = helicopter.b;
            L = helicopter.L;
            f = helicopter.f;
            l_t = helicopter.l_t;
            k1 = helicopter.k1;
            k2 = helicopter.k2;
            lamda_c=helicopter.lamda_c;

            V = mu * omega * R;
            omega_t = V/(mu*R_t);
            A = pi * R^2;
            A_t = pi * R_t^2;

            T = helicopter.thrust(mu)
            C_T = T/(rho*(omega*R)^2*A); % thrust coefficient CT
            sigma = b*chrod/(pi*R);  % rotor solidity
            sigma_t = b*chrod_t/(pi*R_t);  % tail rotor solidity
            C_L = 6*C_T/sigma; % average lift coefficient
            alpha = C_L/a; % average angle of attack
            C_D = 0.0087 - 0.0216*alpha + 0.4*alpha^2; % average drag polar

            D = 0.5 * rho * V^2 * f; % fuselage drag equation
            C_H = sigma * C_D *mu/4; %rotor drag coefficient

            V_1h = sqrt(T/(2*rho*A));
            V1 = V_1h*(sqrt(-0.5*(V/V_1h)^2+sqrt(0.25*(V/V_1h)^4+1))); %hover induced velocity
            lamda1 = V1/(omega*R); % induced velocity inflow factor
            lamda = lamda1 + mu*alpha; % nflow ratio

            C_Q = k1*lamda1*C_T+lamda_c*C_T+sigma*C_D*(1+k2*mu^2)/8+mu^3/2*f/A; % main rotor torque coefficient

            Q = C_Q * rho * pi * R^2 *(omega*R)^2 * R; %torque

            T_t = Q/l_t; % tail rotor thrust coefficien

            C_Tt = T_t/(rho*A_t*(omega_t*R_t)^2); % tail rotor thrust coefficient
            C_Lt = 6*C_Tt/sigma_t;
            alpha_t = C_Lt/a;
            C_Dt = 0.0087 - 0.0216*alpha_t + 0.4*alpha_t^2;

            V_1ht = sqrt(T_t/(2*rho*A));
            V1t = V_1ht*(sqrt(-0.5*(V/V_1ht)^2+sqrt(0.25*(V/V_1ht)^4+1)));
            lamda_t = V1t/(omega_t*R_t); % tail rotor inflow ratio

            C_p = k1*lamda1 * C_T+lamda_c*C_T+sigma*C_D/8*(1+k2*mu^2)+mu^3/2*f/A+(R_t/R)^2*(lamda_t*C_Tt+sigma_t*C_Dt/8*(1+3*mu^2));
            % total helicopter Power coefficient
            y = C_p;
        end
    end
    methods (Static = true)
        % this function is aimed to calculate the required power
        function y = power_required(mu)
            rho = helicopter.rho; 
            R = helicopter.R;
            omega = helicopter.omega;
            A = pi * R^2;

            C_p = helicopter.Cp(mu)
            P = C_p * rho * A *(omega*R)^3;  % watts
            y = P / 1000;
        end
    end
    methods(Static=true)
        % this function is aimed to calculate the climb rate
        function y = climbRate(mu)
            rho = helicopter.rho;
            R = helicopter.R;
            omega = helicopter.omega;
            T = helicopter.thrust(mu);
            Cp = helicopter.Cp(mu)
            k1 = helicopter.k1

            A = pi * R^2;
            V_t = omega * R;
            C_T = T/(rho*V_t^2*A);
            power_a = 1544;
            power_r = helicopter.power_required(mu);
            power_e = power_a-power_r;
            lambda_c = power_e*1000/(rho*A*V_t^3*C_T);
            V_c = lambda_c * omega * R;
            y = V_c
        end
    end
    methods(Static=true)
        % this function is aimed to calculate the climb angle
        function y = theta_climb(mu)
            rho = helicopter.rho;
            R = helicopter.R;
            omega = helicopter.omega;
            T = helicopter.thrust(mu);

            A = pi * R^2;
            V_t = omega * R;
            C_T = T/(rho*V_t^2*A);
            power_a = 1544;
            power_r = helicopter.power_required(mu);
            power_e = power_a-power_r;
            lambda_c = power_e*1000/(rho*A*V_t^3*C_T);
            theta = atan(lambda_c/mu);
            y = theta;
        end
    end
end





